using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics;
using MathNet.Numerics.Transformations;
using System.Drawing;

namespace ParticleTracker.Utilities
{
    public enum PadPattern { DataUpperLeft, DataLowerRight };

    static class XCorr2
    {
        public static Complex[,] xCorr2(Complex[,] image1, Complex[,] image2)
        {
            //This part doesn't work for not the same size pieces.
            //ArrayOps.printArray(image1);
            //Console.Write("\n\n\n");
            //Complex[,] image1z = zeroPad(image1,PadPattern.DataUpperLeft);
            
            //Complex[,] image2z = zeroPad(image2,PadPattern.DataLowerRight);
            //ArrayOps.printArray(image2);
            //Console.Write("\n\n\n");

            int rowdim = nextPow2(2*Math.Max(image1.GetLength(0), image2.GetLength(0)));
            int coldim = nextPow2(2*Math.Max(image1.GetLength(1), image2.GetLength(1)));

            Complex[,] image1z = zeroPad(image1, PadPattern.DataUpperLeft, rowdim, coldim);
            //ArrayOps.printArray(image1z);
            //Console.Write("\n\n\n");
            Complex[,] image2z = zeroPad(image2, PadPattern.DataLowerRight, rowdim, coldim);
            //ArrayOps.printArray(image2z);
            //Console.Write("\n\n\n");

            Complex[,] Fbar = conjugate(fft2(image1z));
            Complex[,] G = fft2(image2z);

            Complex[,] H = new Complex[Fbar.GetLength(0),Fbar.GetLength(1)];

            for (int i = 0; i < Math.Max(Fbar.GetLength(0), G.GetLength(0)); i++)
            {
                for (int j = 0; j < Math.Max(Fbar.GetLength(1), G.GetLength(1)); j++)
                {
                    H[i, j] = Fbar[i, j] * G[i, j];
                }
            }
            Complex[,] h = ifft2(H);

            // trim the array back down to size --
            // get rid of the junk from zero padding
            // for fft.

            // these just happen to be the useful number of items in h
            // in each direction. Draw some pictures to figure out why.
            int rows = image1.GetLength(0) + image2.GetLength(0) - 1;
            int cols = image1.GetLength(1) + image2.GetLength(1) - 1;
            int x = h.GetLength(0) - rows;
            int y = h.GetLength(1) - cols;

            h = ArrayOps.getRegion(h, new Rectangle(x, y, rows, cols));



            return h;


        }

        public static Complex[,] fft2(Complex[,] target)
        {
            ComplexFourierTransformation cft = 
                new ComplexFourierTransformation(TransformationConvention.Matlab);

            Complex[] row;
            for(int i = 0; i < target.GetLength(0); i++)
            {
                //ArrayOps.setRow<Complex>(i, target, cft.TransformForward(ArrayOps.getRow<Complex>(i, target)));
                row = ArrayOps.getRow(i, target);
                cft.TransformForward(row);
                ArrayOps.setRow(i, target, row);

            }

            Complex[] column;
            for (int j = 0; j < target.GetLength(1); j++)
            {
               //ArrayOps.setColumn<Complex>(j, target, cft.TransformForward(ArrayOps.getColumn<Complex>(j, target)));
                column = ArrayOps.getColumn(j, target);
                cft.TransformForward(column);
                ArrayOps.setColumn(j, target, column);
            }

            return target;
        }

        public static Complex[,] ifft2(Complex[,] target)
        {
            ComplexFourierTransformation cft =
                new ComplexFourierTransformation(TransformationConvention.Matlab);

            Complex[] row;
            for (int i = 0; i < target.GetLength(0); i++)
            {
                //ArrayOps.setRow(i, target, cft.TransformBackward(ArrayOps.getRow(i, target)));
                row = ArrayOps.getRow(i, target);
                cft.TransformBackward(row);
                ArrayOps.setRow(i, target, row);
            }

            Complex[] column;
            for (int j = 0; j < target.GetLength(1); j++)
            {
                //ArrayOps.setColumn(j, target, cft.TransformBackward(ArrayOps.getColumn(j, target)));
                column = ArrayOps.getColumn(j, target);
                cft.TransformBackward(column);
                ArrayOps.setColumn(j, target, column);
            }

            return target;
        }

        public static Complex[,] conjugate(Complex[,] target)
        {
            for (int i = 0; i < target.GetLength(0); i++)
            {
                for (int j = 0; j < target.GetLength(1); j++)
                {
                    target[i, j] = target[i, j].Conjugate;
                }
            }
            return target;
        }

        

        public static Complex[,] zeroPad(Complex[,] target, PadPattern pattern)
        {
			int rowDim = nextPow2(2*target.GetLength(0));
            int colDim = nextPow2(2*target.GetLength(1));
			return zeroPad(target, pattern, rowDim, colDim);
			
		}   
            
		
		public static Complex[,] zeroPad(Complex[,] target, PadPattern pattern, int rowDim, int colDim)
		{
			
            //System.Console.Write("rowdim = " + rowDim + "  coldim= " + colDim +"\n"); 
            Complex[,] padded = new Complex[rowDim, colDim];

            if (pattern.Equals(PadPattern.DataUpperLeft))
            {
                for (int i = 0; i < target.GetLength(0); i++)
                {
                    for (int j = 0; j < target.GetLength(1); j++)
                    {
                        padded[i, j] = target[i, j];
                    }
                }
                
            }
            else if(pattern.Equals(PadPattern.DataLowerRight))
            {
                int ii = 0;
                int jj = 0;
                for(int i = rowDim-target.GetLength(0); i < rowDim; i++)
                {
                    jj = 0;
                    for( int j = colDim-target.GetLength(1); j<colDim; j++)
                    {
                        padded[i,j] = target[ii,jj];
                        jj = jj + 1;
                    }
                    ii = ii + 1;
                }
            }
            return padded;
        }
		

        public static int nextPow2(int num)
        {
            int nextPow2 = (int)Math.Pow(2,Math.Ceiling(Math.Log(num, 2)));
            return nextPow2;
        }

        

    }
}
